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Abstract 

We describe the first version (vl.OO) of the code hfbrad which solves the Skyrme- 
Hartree-Fock or Skyrme-Hartree-Fock-Bogolyubov equations in the coordinate rep- 
resentation within the spherical symmetry. A realistic representation of the quasi- 
particle wave functions on the space lattice allows for performing calculations up 
to the particle drip lines. Zero-range density-dependent interactions are used in the 
pairing channel. The pairing energy is calculated by either using a cut-off energy in 
the quasiparticle spectrum or the regularization scheme proposed by A. Bulgac and 
Y. Yu. 
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PROGRAM SUMMARY 



Title of the program: hfbrad (vl.OO) 

Program obtainable from: CPC Program Library, Queen's University of Bel- 
fast, N. Ireland 

Licensing provisions: none 

Computers on which the program has been tested: Pentium-Ill, Pentium-IV 

Operating systems: LINUX, Windows 

Programming language used: FORTRAN-95 

Memory required to execute with typical data: 30 MBytes 

No. of bits in a word: The code is written with a type real correspond- 
ing to 32-bit on any machine. This is achieved using the intrinsic function 
selected_real_kind at the beginning of the code and asking for at least 12 
significant digits. This can be easily modified by asking for more significant 
digits if the architecture of the computer can handle it. 

No. of processors used: 1 

Has the code been vectorised?: No. 

No. of bytes in distributed program, including test data, etc.: 400 kbytes 

No. of lines in distributed program: 5164 (of which 1635 are comments and 
separators) 

Nature of physical problem: For a self-consistent description of nuclear pair 
correlations, both the particle-hole (field) and particle-particle (pairing) chan- 
nels of the nuclear mean field must be treated within the common approach, 
which is the Hartree-Fock-Bogolyubov theory. By expressing these fields in 
spatial coordinates one can obtain the best possible solutions of the problem; 
however, without assuming specific symmetries the numerical task is often too 
difficult. This is not the case when the spherical symmetry is assumed, because 
then the one-dimensional differential equations can be solved very efficiently. 
Although the spherically symmetric solutions are physically meaningful only 
for magic and semi-magic nuclei, the possibility of obtaining them within tens 
of seconds of the CPU makes them a valuable element of studying nuclei across 
the nuclear chart, including those near or at the drip lines. 

Method of solution: The program determines the two-component Hartree- 
Fock-Bogolyubov quasiparticle wave functions on the lattice of equidistant 
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points in the radial coordinate. This is done by solving the eigensystem of two 
second-order differential equations by using the Numerov method. Standard 
iterative procedure is then used to find self- consistent solutions for the nuclear 
product wave functions and densities. 

Restrictions on the complexity of the problem: The main restriction is related 
to the assumed spherical symmetry. 

Typical running time: One Hartree-Fock iteration takes about 0.4 sec for a 
medium mass nucleus, convergence is achieved in about 40 sec. 

Unusual features of the program: none 



LONG WRITE-UP 



1 Introduction 

For several decades, the ground state of nuclei have been studied within the 
self-consistent mean-field approximation. Such a description of the atomic nu- 
cleus, has ability to properly account for the bulk properties of nuclei such as 
masses, radii or shape. The Hartree-Fock method provides a good approxima- 
tion of closed shell magic nuclei, however, the pairing correlations constitute 
an essential ingredient for the description of open shell nuclei. These effects are 
usually described by the Hartree-Fock plus BCS (HFBCS) or Hartree-Fock- 
Bogolyubov (HFB) methods. 

Two main classes of numerical methods have been so far used for the im- 
plementation of the HFBCS or HFB methods. In the first one, one employs 
the expansion of the quasiparticle states on a discrete basis of orthogonal 
functions, usually provided by the harmonic oscillator potential. In this case 
the non-linear HFB(CS) equations are formulated in a matrix form and can 
be solved by using an iterative procedure (hereafter Bogolyubov iterations). 
In the second class, one uses the direct integration of the equations in the 
coordinate representation. This is usually done within the box boundary con- 
ditions to discretize the spectrum of quasiparticle states. The first approach is 
more elegant and efficient for finite-range interactions, however, it has a dis- 
advantage that the use of harmonic oscillator wave functions may perturb the 
correct description of the asymptotics of the system. On the other hand, the 
asymptotic part of the quasiparticle wave functions does not suffer from this 
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possible weakness when the problem is solved in the coordinate basis, but here 
one has to solve a set of integro-differential equations and as a consequence 
this method is usually only applied for the zero-range Skyrme force, for which 
the problem reduces to a set of differential equations. 

The question of the asymptotic properties of the atomic nuclei (neutron skins, 
halos) is an urgent challenge in theoretical nuclear structure, and is particu- 
larly important due to the availability of beams of exotic nuclei ^ , which allow 
for relevant experiments. The HFB problem in coordinate basis, along with 
the use of the effective Skyrme force is indeed a powerful tool for studying 
ground-state properties of nuclei, and such studies have been and will be pur- 
sued especially in the regions of nuclear chart where the experiments are being 
performed. 

Two computer codes that solve self-consistent equations on the basis have 
recently been published [1,2] for non-spherical shapes. We refer the reader 
to these publications for a review of approaches and methods that can be 
and have been used in such cases. Within the assumed spherical symmetry, 
which is the subject of the present study, the first self-consistent solutions were 
obtained in the 1970's [3,4], and these old codes were later used for decades 
by different groups, however, they have never been published. 

The first description and implementation of the HFB problem within the 
coordinate-space spherical symmetry was presented in Ref . [5] . The code writ- 
ten during this work was later updated many times since its early version, and 
distributed widely, however, neither this code was ever the object of an official 
publication. Indeed, the absence of accompanying manual makes it hard to 
use, and its older versions are sometimes used despite the fact that the newer 
ones have been made available. For coordinate-space spherical symmetry, the 
code solving the HFBCS equations for the Skyrme interaction was published 
in Ref. [6] and that for the relativistic mean field Hartree-Bogolyubov method 
in Ref. [7]. A similar unpublished code also exists to treat pairing correlations 
for the finite-range Gogny interaction [8] , although the code that would solve 
the full HFB problem within these conditions is not yet available. 

The aim of this work is to provide a modern code (named HFBRAD) for the 
spherical Skyrme-HFB problem in the coordinate representation. The code 
constitutes a completely rewritten version of that constructed in Ref. [5] , but 
it also features implementation of the pairing renormalization [9]. Although 
such code has to be used with caution in open-shell nuclei, because there the 
deformation effects not included here are essential, it nevertheless provides 
a fast and easy tool for a quick estimate of nuclear properties, which can 
precede much more time consuming calculations beyond the spherical symme- 

^ For example in Europe: ARENAS3, Louvain-la-Neuve; ISOLDE-CERN, Geneva; 
GANIL, Caen; and GSI, Darmstadt. 
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try. In this paper wc also provide a clear description of the different Skyrme 
forces implemented in the code, terms neglected in the energy functional, and 
parameters in the pairing channel. 

This paper is organized as follows. In section 2 we briefly show how the Skyrme- 
HFB equations are derived and discuss their main properties. Notably, a sig- 
nificant part is devoted to the problem of the divergence of the energy due 
to the zero range of the force in the pairing channel. In section 3 we present 
several parametrizations of the Skyrme force implemented in the code, with a 
particular attention paid to the pairing channel. In section 4 we give definitions 
and meanings of the various quantities and observables which are the result 
of calculations. Some aspects of the numerical treatment of the problem are 
presented in section 5, and finally, sections 6 and 7 give a precise description 
of the input and output data files. 



2 The Skyrme Hartree-Fock-Bogolyubov equations 



The HFB approximation is based on the use of a trial variational wave function 
which is assumed to be an independent quasiparticle state |$). This state, 
which mixes different eigenstates of the particle number operator, is a linear 
combination of independent particle states representing various possibilities 
of occupying pairs of single particle states. Following the notations and phase 
convention of [5] we define the particle and pairing density p and p matrices 

by 

p(rag, rVV) = {^al,„,^,a^^q\^) , (1) 



p(r(7g, rVV) = -2a' {^a^>^^,q>a^^q\^) , (2) 

where the operators aj^^ and a^^^q create and annihilate a nucleon at the point 
r having spin a — ±\ and isospin q — :k\. The symmetry properties of p 
and p as well as the relation between p and the pairing tensor k (defined for 
example in [10]) are discussed in [5]. 

The variation of the energy expectation value E = with respect to 

p and p under the constraints N = ($|A^|$) and Z = (<I>|Z|$) (for neutrons 
and protons) leads to the Hartree-Fock-Bogolyubov equation which reads in 
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coordinate representation 



rr' 



h{Ya,v'a') h{ra,r'a') \ (ipi{E,r'a') 
" \ h{ra,r'a') -h{ra,r'a')l \ip2{E,r'a') 




where the particle and pairing fields are given by 

5E 



h{ra,r'a') 



Sp{ra, r'a') ' 



h{ra, r'a') 



6E 



Sp{ra, r'a') 



(3) 



(4) 



Once again we refer the reader to the article [5] for the discussion concern- 
ing the quasiparticle spectrum, its symmetries and the relations between the 
components of the HFB spinors and the densities. 



2.1 Local densities 



In the Skyrmc-HFB formalism, the evaluation of the expectation value of 
the energy leads to an expression which is a functional of the local densities, 
namely, the particle (normal) and pairing (abnormal) densities 

p(r) = "P^i^i^ P(0 = - I] "PiiEi, r)(p2{Ei, r) , (5) 

i i 

and their derivatives. The presence of non local terms in the force leads to a 
dependence on the normal and abnormal kinetic densities 



r{r)^j:\VMEi:r)\\ (6) 

i 

f(r) = - E V(^i {E„ r) • Vip2{Ei, r) , (7) 

i 

while the spin-orbit term leads to a dependence on the spin current tensors 
and Jij . We do not give here the definitions of these tensors, because for 
the spherical symmetry discussed here they reduce to the corresponding spin 
current vectors J and J 

J(r) = i J2 ¥>2{Ei, v)V^2{Ei, r) {a'\&\a), (8) 

i 

J(r) = -i ^ (^i(^,, v)V^2{Ei, v){a'\&\a). (9) 
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2.2 Spherical symmetry 



In the code HFBRAD the solutions are restricted to have a spherical symmetry 
and do not mix proton and neutron states. In this situation the wave functions 
have the good quantum numbers {nijmq) {n is a good quantum number since 

the spectrum is discretized inside a spherical box) and all the solutions inside 
an (ni?jg)-block are degenerated. Furthermore, the radial part of the wave 
functions can be chosen to be real. Thus we use the ansatz 



^E,ra) = !lit!lAZlYW(f)(£m4abm) , t = 1, 2 , (10) 



for the wave functions. The local densities can be written using the radial 
functions (omitting the isospin quantum number) 



p(r) 



p(r) 



^(2j + l)«2(n£j,r), 



nij 



47rr^ 



^(2j + l)ui{nij,r)u2{nej,r) . 



nij 



;i2) 



For the kinetic densities we have 



r r 



E 



nij 



2j + l 



' U2{nij,r) Y , £{£+!) , 
u^{nij, r) H — u^irdj, r 



2j + l 
47rr^ 



+ 



nij 



u[{nij,r) - 



ui{nij,r) 



u'^inij^r) - 



Ui{nij,r)u2{nij,r) 



, (13) 
U2{n£j,r) '^ 

(14) 



Finally the spin current vector densities have only one non vanishing compo- 
nent given by 



J(r) 



J(r) 



E(2J + 1) 



nij 



E(2J + 1) 



nij 



j{j + I)- I) - 



ul{nej,r), (15) 
ui{n£j, r)u2{n£j, r|16) 
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2.3 The Hartree-Fock-Bogolyubov energy 

In the Skyrme-HFB approximation, the total energy £^ of a nucleus is given 
as a sum of kinetic, Skyrme, pairing and Coulomb terms: 

E — K + -Bskyrme + -E'pair + EqouI 

=y dh [/C(r) + £skyrme(r) + ^pair(r) + ^Coul(r)] • (17) 

The derivation of the energy is explained in detail in several articles (see Refs. 
[3,11] for the HF energy and [5] for the HFB case). The expression of the 
Skyrme force is given in section 3. 

In equation (17), the kinetic energy of both neutrons and protons is given by 
isoscalar kinetic density. The neutron and proton masses being approximated 
by their average value, and one has 

where the factor in parentheses takes into account the direct part of the center- 
of-mass correction [12]. Since wc only consider cvcn-cvcn nuclei, the Skyrme 
part of the energy functional is time even, it can be written as a sum of 
isoscalar (T = 0) and isovcctor (T = 1) parts (see e.g. Ref. [13]) or as a sum 
on isoscalar, neutrons and protons densities. 



^Skyrme — 



+ 



tl 



1 + 
1 + 



Xi 



^0 + 2 IIIa 



^2 + ^)E 



Pq^q - \ {"^Pqf 



(19) 



Index q stands for neutrons and protons while the absence of index indicates 
the total (isoscalar) density. Using the same notations for the densities, the 
pairing energy density reads 
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•pair — 



(20) 



In this last expression we have added a prime to the parameters of the interac- 
tion. Although the derivation of the general Skyrme energy density functional 
is based on a unique force, its effective nature justifies the use of different sets 
of parameters in the particle-hole and particle-particle channels. The prime 
indices anticipates this possibility which is discussed in section 3. 

The energy density defined in (17) involves a Coulomb term for protons. This 
term contains a direct part which can be expressed using the charge density 
and leads to a local field after variation of the energy and an exchange part 
which would lead to a non local potential. Both parts deserve a special dis- 
cussion. The direct part of the Coulomb energy depends on the charge density 
Pc?i(i') and reads 



We use the point proton density to simplify this expression, i.e. the charge 
density is replaced by the proton density Pp(r). Nevertheless, the proton form 
factor is taken into account when we calculate the charge radius of the nucleus 
(see eq. (57)). The exchange part of the Coulomb energy leads to a non local 
term and is treated with the Slater approximation. This approximation means 
that we keep only the first term of the density matrix expansion in the local 
density approximation [14] 



The error introduced by this approximation has been estimated by study- 
ing the next order term [15] or more recently by comparison with the exact 
treatment of the Coulomb energy [16]. The Slater approximation seems to 
have little consequence on the nuclei bulk properties although some signifi- 
cant effects can be expected on the position of the proton drip-line and the 
Coulomb displacement energy of single particle levels. Finally, the Coulomb 
contributions to the pairing energy and fields are not included in the HFBRAD 
code. 




(21) 




(22) 
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2-4 The Hartree-Fock-Bogolyubov mean fields 



Since the Bogolyubov transformation does not preserve the particle number, 
we introduce two Lagrange multipUcrs Ajv and \z to conserve the average 
neutron and proton number. The HFB equations are then obtained by writing 
the stationary condition 6[S — {\nN + \zZ)] = 0. The dependence of S on 
the kinetic densities leads the effective mass 



2m* 

2m 4 



1 + 



Xi 



l + f)p+(x, + i)p. 



and to the abnormal effective mass 



Mr, 



|(l-^'l)P. 



(23) 



(24) 



(25) 



The particle-hole (Hartree-Fock) fields is given by 



f)p-(^o + ^)P, 



+ 



1 + 



Xi 



-I 

t' 



2 

,T3 



2 ' J 

(2 + 7)p^+^ 



2^3 + 



1 



IP 



,7-1 



9' 



24 

-^(VJ + VJ,) 

In the case of protons, varying the expressions (21) and (22) leads 
following expression for the Coulomb field 



K(r) = ^ 



'^-e^(^)%fW. 



(26) 
to the 

(27) 



The particle-particle (pairing) field is 

Uq = 1(1 - A)Pq + t(i - ^'i) h - Wp<^ 



|(l-4)p^'p,. (28) 
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Finally, the spin-orbit fields (oc ^ • s) have the following form factors: 



^i^'l (^1^1 + ^22^2) + ^(^1 - t2)Jg + WoW{p + pg) , (29) 



|(i + 4) + w^ 



Jn 



(30) 



2.5 The Hartree-Fock-Bogolyubov equations 
Writing the fields in matrix form 
M = 




The HFB equations read 



ar ar r 







= E 










w 



(31) 



(32) 



(33) 



An r-dependant mixing of components and scaling described in Ref. [5] allows 
us to write this equation as an equation with no differential operator in the 
coupling terms and no first order derivative 



M*i,^ + Vh + Wh = Eh, 



dr 

M*£,h-Vf2 + Wh^Eh. 



(34) 



This last form with no first order derivative of the functions is particularly 
suitable for the numerical integration by the Numerov algorithm briefly dis- 
cussed in section 5.2. 



2.6 Asymptotic properties of the HFB wave functions 



The asymptotic properties of the two components of the HFB quasiparti- 
cle wave functions, and their dependence on A and E, were discussed in 
Refs. [17,5,18]. Here we complement this discussion by further elements, which 
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pertain mainly to weakly bound systems where the Fermi energy A is very 
small. 

Assuming that M*—-^, which anyhow is always fulfilled in the asymptotic re- 
gion, and neglecting for simplicity this trivial mass factor altogether, Eqs. (34) 
have for large r the following form 

/f +(E + A)A= Wh, , , 

(35) 

Here we have also neglected all mean-field potentials, assuming that the dis- 
tance is large enough that they are much smaller than the two characteristic 

constants k"^ — E + X and — E — \. Moreover, wc consider the case of small 
Fermi energy such that there are no discrete HFB states, i.e., E>—X. Within 
such conditions, the only remaining question is whether in the asymptotic 
region we can neglect the couphng potential W. 

In order to discuss this question, we note that the coupling terms can be con- 
sidered as inhomogeneities of the linear equations (35), and therefore, asymp- 
totic solutions have the form 

/i(r) oc sin(A;r + 5) + Fi(r), 

(36) 

/2(r)(x e— + F2(r). 

Wc can now iterate these solutions starting with Fi{r)=F2{r)=0, which cor- 
responds to neglecting the coupling terms in the zero order. Then, in the first 
order we see that in the equation for fi the coupling term can still be ne- 
glected as compared to the large term {E + X)fi. Therefore, to all orders we 
have Fi(r)=0, and the asymptotic solution /i(r) oc sin(A;r -|- 6) defines the 
inhomogeneity Wfi of the equation for /2. 

We arrive here at the conclusion that the asymptotic form of may involve 
two terms and a more detailed analysis is needed before concluding which one 
dominates. To this end, we note that the coupling potential W depends on the 
sum of products of lower and upper components, and therefore has a general 
form of 

W{r) oc W^osc(r)e-^7r2, (37) 

where Wosd^) is an oscillating function of order 1, is the decay constant, 
which may characterize another quasiparticle state than the one with energy 
E, and the factor comes from the volume element, cf. Eq. (10). We see here 
that the term F2{r) in the asymptotic form of f2{r) vanishes as e"^*", but it is 
difficult to say which one of the two terms dominates asymptotically. We only 
note that for small quasiparticle energies the decay constant k, is also small 
and hence the first term is then more likely to dominate. 
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5 10 15 20 25 30 5 10 15 20 25 30 

r [fm] r [fm] 



Fig. 1. HFB quasiparticle wave functions ui{r) and U2{r) of the four lowest states 
in i^^Sn. 



In order to illustrate the above discussion, we have performed the HFB calcu- 
lations in -'^'^'^Sn, where A„=— 0.345 MeV and the lowest £=0 quasiparticle state 
of i5=0.429MeV (for i?box=30fm) leads to a very diffused coupling potential 
with small decay constant Quasiparticle wave functions corresponding to 
the four lowest £=0 quasiparticle states are shown in Fig. 1. One can see that 
the asymptotic forms of the second components U2{r) of the two lowest quasi- 
particle states are not affected by the second term F2{r), at least up to 30 fm. 
Only the third and fourth states switch at large distances to the oscillating 
asymptotic forms with a smaller decay constant. This happens at rather large 
distances where the densities are anyhow very small. Hence the change in the 
asymptotic properties does not affect any important nuclear observables. 

In all cases that we have studied, the practical importance of the second term 
F2(r) is negligible. However, its presence precludes simple analytic continua- 
tion of the wave functions in the asymptotic region. We would like to stress 
that the asymptotic forms discussed in this section are numerically stable, 
and that they are unrelated to the numerical instabilities discussed in Sec. 5.3 
below. 
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2. 7 Pairing correlations and divergence of the energy 

The use of a zero range interaction in the particle-particle channel leads to 
the energy divergence problem [9,19] related to the fact that the short-range 
behaviour of the abnormal density is 

pir) = lim p(r,r') oc -j -. (38) 

r'-^r r — r 

The other abnormal densities diverge as well, but along that the pairing cor- 
relations also induce the divergence of the kinetic density r(r). One simple 
way to overcome the divergence of the energy is to introduce a cut-off energy 
E'cut in the summations of all the densities 

E< - E< • (39) 

Ei>0 0<Ei<Ecut 

This cut-off somehow simulates the finite range of the interaction in the pairing 
channel. It therefore constitutes an additional parameter of the interaction 
and must be adjusted along with the other quantities which define the pairing 
interaction, see e.g. Ref . [18] . 

A more elegant way to prevent the divergence was recently proposed in Refs. 
[9,19]. It is based on the subtractions of the divergent part of the abnormal 
density. In its present implementation in the code HFBRAD, this method can 
only be applied if the pairing force does not contains non-local terms {i.e. 
gradient terms). In that case, the pairing energy depends only on p (and 
possibly on p but not on f or J), i.e., 

pI = T.9[p]pI (40) 
q 

and the pairing potential depends linearly on p, 

U,^2g[p]p,. (41) 
The regularized pairing field can then be written as 

Ul^^ ^ 2g[p]pl^^ ^ 2gMP, ■ (42) 

Calculation of the various densities is performed by summing up the quasi- 
particle contributions up to a maximum energy E'max, but because of the reg- 
ularization, this maximum energy has not the meaning of a cut-off energy 

Ecut- 

It can be shown [19] that when we evaluate the total energy of the system. 



pair 



E 
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the pairing energy must be defined as 

^pair = l]c/efr[p]Pg- (43) 

With this definition, the kinetic and pairing energies both diverge but the 
sum, and indeed the total energy, does not [19]. If the pairing energy density 
involves kinetic or gradient terms (like it was implemented for the Skyrme 
interaction SkP [5]), the regularization would require the subtraction of higher 
order terms. This possibility in not implemented in the present version of 
the code. Finally, let us note that relations (21)-(26) in Ref. [19] imply that 
the variation of quantities kp and kc with respect to the densities should 
be included in the mean fields. This contribution is supposed to be small 
and is not taken into account in the code. On the other hand, for density- 
dependent pairing interaction, the variation with respect to the explicit density 
dependence is taken into account in the code, and gives a contribution to the 
energy of the order of a few keV. 



3 The effective Skyrme force 

For a general overview of the foundations and properties of the Skyrme force 
we refer the reader to the review article [20] and references therein. The Skyrme 
force is an effective interaction depending on a limited number of parameters, 

Vi2 ^to{l + XoP^)5 + + XiP^) (k''5 + 5k^) + ^2(1 + X2Pa)\i' ' Sli 

+ + x^P,)p^5 + iVFo (o-i + a2) ■ (k' X 5k) , 

where 5 is a short notation for 5{yi — r2), k = ki — k2 acting on the right and 
k' = ki — k2 acting on the left. 

The parameters of the Skyrme forces were fitted in the literature to reproduce 
various bulk nuclear properties as well as selected properties of some nuclei 
(usually doubly magic nuclei). Simplifications have often been made in the 
expression of the functional (19), like treatment of the Coulomb exchange term 
with the Slater approximation and/or omission of the two-body center of mass 
contribution or of the " " terms. The latter corresponds to the fourth line 
in Eq. (19) and to the two first terms in the spin-orbit mean field, Eq. (29). It 
is important to keep in mind that one should use each parametrization of the 
functional within the same simplifications with which it has been adjusted to 
data. 

For all forces implemented in HFBRAD, the force in the particle-particle channel 
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is chosen to be 

W2=(t'o + |p^')5- (45) 

Table 1 gives the different sets of parameters for the force in the particle-hole 
channel, while table 2 gives the parameters in the particle-particle channel. It 
is important to keep in mind that these latter parameters have been adjusted 
along with a given cut-off and cut-off diffuseness, which are two additional 
parameters of the force. 





SIII [4] 


SkM* [21] 


SLy4 [22] 


SLy5 [22] 


SkP [5] 


to 


-1128.75 


-2645.0 


-2488.913 


-2488.345 


-2931.6960 


h 


395.0 


410.0 


486.818 


484.230 


320.6182 


t2 


-95.0 


-135.0 


-546.395 


-556.690 


-337.4091 


t3 


14000.0 


15595.0 


13777.0 


13757.0 


18708.96 


Xq 


0.45 


0.09 


0.834 


0.776 


0.2921515 


Xl 


0.0 


0.0 


-0.344 


-0.317 


0.6531765 


X2 


0.0 


0.0 


-1.000 


-1.000 


-0.5373230 


X3 


1.0 


0.0 


1.354 


1.263 


0.1810269 


7 


1.0 


1/6 


1/6 


1/6 


1/6 


Wo 


130.0 


120.0 


123.0 


125.00 


100.000 


J2 


No 


No 


No 


Yes 


Yes 



Table 1 

Parameters in the particle-hole channel for the different versions of the Skyrme 
forces implemented in the code HFBRAD. The last line indicates if the terms are 
included in the Skyrme functional when the force is used. 



Three kinds of pairing forces have been adjusted: (i) the volume pairing (the 
pairing field follows the shape of the density), below denoted by "p", (ii) the 
surface pairing (the pairing field is peaked at the surface and follows roughly 
the variations of the density), below denoted by "5p", and (iii) the mixed 
pairing (a compromise between the two first two), below denoted by "p+5p". 
All pairing parameters have been adjusted in order to give a mean neutron 
gap of 1.245 MeV in ^^ogn. 
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cut-ofF 


rej 


jularization 




P 


5p 


p + 5p 


V 

P 




{p + SpY 


t()(SIII) 


-159.6 


-483.2 


-248.5 


-197.6 


-807.0 


-316.9 


t'o(SkM*) 


-148.6 


-452.6 


-233.9 


-184.7 


-798.4 


-300.7 


t^(SLy4) 


-186.5 


-509.6 


-283.3 


-233.0 


-914.2 


-370.2 


io(SLy5) 


-179.9 


-504.9 


-275.8 


-222.7 


-901.9 


-356.6 


t'o(SkP) 


-131.6 


-429.5 


-213.1 


-196.6 


-1023.0 


-326.5 







-37.5 


-18.75 1() 





-37.5 1[) 


-18.75 1[) 


y 


1 


1 


1 


1 


1 


1 



Table 2 

Parameters in the particle-particle channel for the different versions of the Skyrme 
forces implemented in the code HFBRAD. The three first columns corresponds to the 

use of a cut-off at 60 MeV with a Fermi shape and a diffuscncss of 1 McV, the three 
columns on the right corresponds to the regularization procedure of the pairing field 
described in Sec. 2.7. 



4 Observables and single particle properties 



^.1 Hartree-Fock equivalent energies, radii, and nodes of quasipartide wave 
functions 



For each quasipartide state, in addition to its quasipartide energy Ei given 
by the HFB equation, one has several other characteristics, which are defined 
the following way: 

• the occupation factor Ni which is the norm of the lower component of the 
HFB quasipartide wave function; 

• the (Hartree-Fock) equivalent energy Si which is defined by applying the 
ECS type formula to the occupation factor N^, i.e.. 




which gives 

£, = A + ^,(2iV,-l) ; (47) 

• in the same spirit, the equivalent gap A, is defined by applying the BCS 
type formula to the HFB quasipartide energy Ei, i.e., 

Ei = ^{e, - A)2 + A? , (48) 
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which gives 

A, = 2Ei^Ni{l - Ni) ; (49) 

• the rms radius ^JJrf) is calculated using the lower component of the quasi- 
particle wave function 

(r^)- I rV|,(r)ci^r; (50) 

• the number of nodes of the HFB quasiparticle wave function is defined as 
the number of nodes (including the node at the origin but not the one at 
infinity) of the component which has the greatest amplitude. 

These characteristics of the quasiparticle states can be found in the output 
file hfb_n_p.spe, see Sec. 7. 



4-2 Canonical basis 



The canonical states are obtained by diagonalizing the density matrix (see e.g. 
Refs. [10,18] for the interpretation of the canonical basis) 

jdh'p{vy),p,{v')=vl,Pi{v). (51) 

As discussed in Refs. [5,18], all canonical states have localized wave functions 
and form a basis. The energies of the canonical states are defined as the 
diagonal matrix elements of the Hartree-Fock field h in the canonical basis, 

Si={A\h\A). (52) 

and the pairing gaps associated with these states are the diagonal matrix 
elements of the pairing field, 

= {i^^^h\i^^) . (53) 

Finally, the canonical quasiparticle energy can be defined using the same kind 
of formula as (48) , 

Er = , (54) 

while the occupation probabilities of canonical states are given by 

1 / _ 

These characteristics of the canonical states can be found in the output file 
hfb_n_p.spe, see Sec. 7. 
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4-3 Observables and other characteristic quantities of the system 



The rms radii for protons and neutrons are defined as {q — n or p) 



) = j^ r'pl{r)dh. (56) 



The charge radius is obtained from the proton radius by taking into account 
the proton charge distribution in an approximate way, i.e., 

K) = {rl) + {r)l. (57) 
with (r) p — 0.8 fm. The mean gaps are the average values of the pairing fields 

(A,> = -^^. (58) 

The fluctuations of the particle numbers are defined as {N^ — {Nq)'^) and given 
by 2Tr [pI - Pq . 

Finally, the rearrangement energy, which comes from the density dependence 
of the force, and which shows how much the force is modified by the medium 
effects, is given by 



TTI ^3 ^ 



^3 



+|v(i-4)/i:p?-?(^)%f 



(59) 



5 Numerical treatment of the problem 



5.1 The problem inside a box 



The system of differential equations (34) is solved in coordinate representation 
in a spherical box for a given choice of boundary conditions. Here we briefly 
introduce notations that useful for the discussion that follows and discuss 
the validity of approximations to solve the differential problem. The radial 
coordinate r is approximated by a mesh of points spaced by the step h: r„ = 
nh, with n G {0, 1, ...,N}. For any quantity F(r) approximated on the mesh 
we use notation F„ = F{r-n)- 
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On the wall of the box, wc use two possible boundary conditions: vanish- 
ing wave functions (Dirichlet condition) or vanishing derivatives of the wave 
functions (Neumann condition). With this choice, the continuous part of the 
spectrum is discretized. Important questions are how well the continuum is ap- 
proximated by this discretization and whether it does not introduce numerical 
artifacts when solving system of equations (34). A discussion concerning the 
behaviour of the resonant like solution inside the box can be found in [5]. The 
non-linearity of the equations makes it difficult to discuss further the effect 
of the discretization of the spectrum, nevertheless, the numerical examples 
shown in Sec. 8 show that such numerical method provides a fairly good ap- 
proximation of the exact problem. In particular, it can be observed that once 
a sufficiently large box is used, the global properties of the nucleus become 
perfectly independent of its size. 



5.2 The Numerov algorithm for the HFB equation 



The Numerov (Cowell) method [23] is the standard technique for a precise in- 
tegration of second-order differential equations; here we only give some details 
pertaining to its application to the system of two HFB equations (34). Let 
function y(r) be the solution of the differential equation 

/ = F(y,r). (60) 

The Numerov algorithm is based on the finite difference formula for three 
consecutive points on a mesh with step /i, 

- 22/„ + = ^ + lo^:: + y':,^^) + oih'^), (ei) 

where o{h^) means that the error in this relation is of the order h^. Combining 
this formula with (60), one obtains the relation 

( 1 - ^^n+l j 2/n+l - 2 h + —FA Vn+il- ^^i^n-l j 2/n-l = 0, (62) 

which is obviously also valid if y is a vector and F a square matrix. Introducing 
notations 

= (l - , (63) 

and 

Gn = AnVn , (64) 
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one obtains the recurrence relation (Numerov iterations), 



Gn+l — 12yn ~ lOGn ~ G'n-l , 



(65) 



In the present apphcations, A„ is a 2 x 2 matrix (position and energy depen- 
dent) and the calculation of its inverse could be the bottleneck of the procedure 
if it is not computed cleverly. This can be done in the following way. Using 
the notations introduced in Sec. 2.5, one can write matrix An in an explicit 
form as 



I _ h^Va_ I h'^ E 



12 M* 



h?_Wn 
12 M* 



\ 



\ 



12 M* 



whereby its inverse is given by 



12 M* VI Ml) 



(66) 



1- 



12M* 



fe2 
12M* 



V 



1 _ l^Yn. _ E 

12 M* 12 M* 



12 M* 



h^Wn 1 _ hiVn_ I E 

12 M* ^ 12 M* 12 M*J 



It appears clearly that a lot of time can be saved by calculating and storing 
the energy independent terms. 



12M* ' 12M* ' 12M* 



and 1 



12M* 



+ 



JfWn 

12M* 



(67) 



only once in each (£jg)-block; then the energy-dependent terms can be in each 
iteration calculated very rapidly. 



The eigen energies are found by integrating two linearly independent solutions, 
and y*^^+)(r), from the origin to a given point r^, called the matching 



point. These two solutions are selected by imposing that either the first or the 
second component of the quasiparticle wave function vanishes near the origin, 
i.e., 

/ 



whereby we have 








for r — > , 



(2+) 



(68) 



(69) 
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and 



(1+) 




yi 



(2+) 




(70) 



Next, another two linearly independent solutions, y 



(1- 



and y 



are 

found by the backward integration from the wall of the box i?box = f^N r^^, 
again imposing that either the first or the second component of the quasipar- 
ticle wave function vanishes near for the Dirichlet boundary conditions, 
i.e.. 



(1-) 



y 



whereby we have 



and 



V 







y 



(2-) 



for r —>■ , 



J 



yV 



^(1-) _ 

i/Ar-i — 




(2-) 

yN 



(2-) 




(71) 



(72) 



(73) 



For the Neumann boundary conditions, at the left-hand-sides of Eqs. (71)- 
(73) one replaces functions by derivatives. Finally, one obtains four solutions 
valid everywhere except at the matching point and depending on the four 
constants A, B, C, and D. 

The boundary conditions (69)-(70) and (72)-(73) determine the initial values 
(Go, Gi) and (Gat, Gn-i)-, respectively, needed to start the Numerov iterations 
of (65). Obviously, a correct limit should be taken when calculating values of 
Go from Eq. (64). It turns out that one obtains Go=0 for all partial waves 
except I — 1, which is the case that deserves a little more attention. For 
example, focusing our attention on one of the solutions. 




(74) 



we use the fact that in the vicinity of the origin A{r-) is dominated by the 
centrifugal term, so that its diagonal matrix elements are close to This 
leads to the result 



gJ^+^ = lim A(r)y(r) = 



AH" 1 1 



(75) 



An analogous result can be obtained for the other solution at the origin. 
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The eigenenergies are found by requiring that the full solution y — Ai/i^^ + 



-Cy\ 



+ Dy2 ^ and its first derivative are continuous at r = (see 



for example [24]). This matching conditions read 



/,,(!+) ,,(2+) 

ilm tlm 

. /(!+) /(2+) 

\y m y m 



y m 



_y(2-) 
y m 



(a) 






B 




B 




= M4 




C 




C 









0. 



(76) 



The HFB quasiparticle wave functions having two components, the matching 
conditions amount to searching for zeros of the determinant of the 4x4 matrix 
M4 in function of the quasiparticle energy. For £ — 0, the matching point 
is chosen to be near the point where the central part of the Hartree-Fock 
potential is half of its maximum, while for increasing values of £ it is gradually 
shifted to the exterior. A systematic search allows to find the intervals where 
the zeros of this determinant are located, and then the Newton method is used 
to find their exact locations to an arbitrary precision. 



For vanishing determinant, Eq. (76) is used to determine three of the multi- 
plicative constants in terms of the fourth one. This is done by extracting a 
3x3 matrix M3 out of M4 and inverting it. Among the four possible choices 
one may have for M3 we retain the one for which the product M3 x M^^ 
is closest to unity. Finally, the fourth multiplicative constant is fixed by the 
normalization condition of the wave function. 



When approaching the self-consistent solution, the eigenenergies do not change 
very much between two consecutive Bogolyubov iterations. We take advantage 
of this fact by keeping in memory the intervals where the solutions are located 
and using them as initial guesses for the next iteration. This results in finding 
solutions for eigenenergies faster and faster as we get closer to the converged 
self-consistent state. 



5.3 Numerical instabilities 



If we try to solve the HFB equations in rather big boxes (typically beyond 30 
or 40 fm depending on the mass of the nucleus) we encounter strong numerical 
instabilities. This problem is not trivially related to the accumulation of errors 
during the Numerov iteration, because it disappears if the nucleus is treated 
within the HF approximation instead of HFB. The origin of the problem is in 
the fact that at large distances the upper and lower HFB wave functions differ 
by many orders of magnitude, and thus calculations within any number of 
significant digits must fail at some distance, especially if the small component 
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decreases very fast exponentially. This is so because the small component 
becomes easily polluted by the large component even if they are coupled by 
very small matrix elements. 

This numerical instability can be removed by neglecting the pairing field be- 
yond a certain large distance i?cut- Instead of the full pairing field we consider 
the truncated field given by 



Acut(r) = 



It will be shown in section 8.1 that this truncation has no significant effects 
on the solution if i?cut is far enough from the origin. Note that by removing 
the coupling between the lower and upper components one also removes the 
second term in the asymptotic form of the lower component, see Eq. (36). 

Another source of numerical instabilities may come from an accidental degen- 
eracy of deep hole and particle states. In the HFB approach the deep hole 
states are embedded in the continuum of particle states. The deep hole states 
are extremely narrow and thus insensitive to the size of the integration box. 
On the other hand, the discretized particle states are driven by the boundary 
condition and their spectrum becomes denser as the box is enlarged. As a 
result, a deep hole state and a quasiparticle (particle like) state can become 
quasi degenerate and the systematic search for the solution [corresponding to 
the matching condition (76)] can miss both of them simultaneously. This situ- 
ation can show up during the iterations and disappear at the final solution; in 
such a case this numerical accident has no consequence on the final solution. 
If the quasi degeneracy is met close to the final self-consistent solution, the 
iterations may not converge. This can be overcome by using a smaller energy 
step or by slightly changing the size of the integration box in order to lift the 
degeneracy. 



A(r) for r < i?cut , 
for r > i?cut ■ 



(77) 



6 Input data file 



The code HFBRAD is driven using a single input file which contains all infor- 
mation on the geometry of the system, parameters of the force, demanded 
output and nuclei to be computed. The data are read in namelist, which is a 
Fortran-90/95 feature. All variables in a namelist are optional and acquire a 
default value when absent. 
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6.1 The namelist input 



The namelist input contains the variables summarized on table 3 with their 
default values. The variables are 

• force: sets the parameters of the force in the particle-hole channel, pos- 
sible values are SIII, SKM*, SKP, SLY4 and SLY5, it must be written in 

uppercase. 

• mesh_points: number of mesh points in the box minus one since integration 
is made from i — to i — mesh_points. 

• integ_step: distance between two points on the mesh, in fermis. 

• itmax: maximum number of iterations to obtain the convergence of the 
mean fields and total energy. 

• bogolyubov(2) : logical variables to choose between the HF or HFB approx- 
imation, the first element is for neutrons and the second one for protons. 

• eps_energy: maximum relative variation on the total energy before the 
iterations stops. 

• max_delta: in addition to the previous convergence condition, the maximum 
variation on the sum of the neutron and proton pairing gaps must be smaller 
than this value. 

• regularization: logical variable used to switch on or off the regularization 

of the abnormal density. 

• pairing_f orce: integer used to choose the pairing interaction between: full 
Skyrme force (0), volume pairing (1), surface pairing (2) or mixed pairing 
(3). 

• boundary_condition: integer, Dirichlet condition at the box radius (0), 
Neumann (1), Dirichlet for even i and Neumann for odd £ (2), or Neumann 
for even £ and Dirichlet for odd i (3). 

• xmu: part of the field of the previous iteration which is kept to build the 
field for the next iteration. 



6.2 The namelist miclens 



This namelist describes the nucleus to be computed, it can be repeated as 
many times as needed in the input file. Along with the neutron and proton 
numbers, other quantities can be set, notably the parameters of the Skyrme 
force in the particle-particle channel. If a variable is not present in the list and 
has not been present in the previous namelists nucleus, it acquires a default 
value (see table 4). If a variable is not present but was present in a previous 
namelist, it keeps its previous value. The variables are 
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V Ctl IdL/iC 


lype 


X-'ClClLlit VCtlLlC 


force 


chara.cter*4 


GT V/I 


inesli_poiiit s 


integer 


oU 




rea<l 


n 9R 
u . zo 


itiniax 


integer 


boU 






C /T T/1 


Q T\ C Q Q T" m T 

CPb_CllCX 


reBil 




iii8LX_d.e Xt a. 


real 


o . e 1 


regularization 


logical 


r 


pairing_f orce 


integer 


1 


boundary_condition 


integer 





xmu 


real 


0.8 



Table 3 

Variables from the namelist input. The type real means real of the kind chosen in 
the module in cste.f90. 

• neutron and proton: number of neutrons and protons for the nucleus. If 
one of this number is negative then the program will let vary the number 
of particles in order to reach the corresponding drip line. 

• j_max(2): maximum values of 2 J (first for neutrons, second for protons). It 
is reasonable to use different values for neutrons and protons in the case of 
very asymmetric or semi magic nuclei. 

• cut_of f and cut_dif f useness: If the regularization procedure is used, the 
densities are built by summing up contributions from the quasiparticle states 
up to the maximum equivalent energy (47) £^max given by cut_off , if not, 
the summation is made up to the cut-off equivalent energy E'cut (47) given 
by cut_of f with a smooth Fermi profile of diffuseness cut_dif fuseness, 
see Sec. 2.7. 

• r_cut: value of distance i?cut beyond which the paring fields are neglected 
when integrating the HFB equations, see Eq. (77). 

• e_step: initial step in energy (in MeV) for the search of the solution of the 
HFB equations. 

• sktOp and sktSp: parameters of the force in the particle-particle channel 
as defined in equation (45). 

• read_pot: This string can be used to give a name of a file where the initial 
potentials are read. The number of mesh points can be different in the run 
and in the saved potentials but the integration step must be the same. 

• densities, meanfields, quasiparticle and canonical_states: logical 
fiags which can be set if one wants to save the densities, mean fields and 
quasiparticle and canonical states wave functions at the end of the run. 
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VdilidiUit; 


type 




neutron 


integer 




proton 


integer 




j_max(2) 


integer 


(/21,21/) 


cut_off 


real 


OU.U 


cut_diffuseness 


real 


1 n 
i.U 


r_cut 


real 


oU.U 


e_step 


real 


calculated 


sktUp 


real 


predefined 


sktop 


real 


preaerinea 


ifc!cxQ_pOL 


/^T^ o T' o /~i"v ^ /I 




densities 


logical 


XT' 

r 


meanfields 


logical 


F 


quasiparticles 


logical 


F 


cauoiiicaLslaU-s 


loojcal 


F 



Table 4 

Variables from the namelist nucleus. The neutron and proton numbers have no 
default values. If the energy step is not present in the list it is calculated as a 
function of the box radius. By default, the values of sktOp and sktSp are fixed by 
the kind of pairing force chosen in the namelist input. 



7 Output files 



In addition to the optional files containing the densities, mean fields and wave 
functions, at the end of the run the following output files are saved: 



• hfb_n_p.spe: where n and p are the neutron and proton numbers. This file 
contains the quasiparticle spectrum, the canonical states spectrum (if the 
flag canonical_states was set) and the mean fields. 

• hf b . summary: this files summarizes the global properties for all the nuclei 
in the run. 
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Boundary condition 


E (A^) \m ^{r%) 


Dirichlet 
Neumann 


-1131.863 -1.431 -1.067 5.264 
-1131.862 -1.431 -1.067 5.264 



Table 5 

Total binding energy, mean pairing gap, Fermi energy (all in MeV) and neutron 

rms radius (in fm) in ^^'^Sn, calculated for the force SLy4''^^'' and two different 
boundary conditions: Dirichlet (vanishing wave functions at r = i2box) or Neumann 
(vanishing derivative of the wave functions at r = -Rbox)- 

8 Examples 

8.1 Self- consistent calculation for ^^^Sn 

As a first example we discuss here the results obtained for the semi magic 
isotope i^°Sn {N = 100 and Z = 50) using the effective force SLy4''+''''. This 
nucleus is an interesting example since it has a rather small neutron Fermi 
energy (Aat ~ —1.066 MeV). Small neutron separation energy allows us to 
discuss possible effects of the integration box size as well as the choice of the 
boundary condition. If not explicitly stated otherwise, in all calculations the 
integration step is h — 0.2 fm for the box size of i?box = 30 fm, i.e., the 
number of integration points is 150+1. The global properties of the nucleus 
for Dirichlet or Neumann boundary condition are shown in table 5. 

In addition to the small neutron separation energy, one finds several single 
particle states with a significant occupation around Xjy with different values 
of i and indeed different asymptotic bahaviour. The single particle neutron 
states in the vicinity of the Fermi energy are shown in table 6. It is gratifying 
to observe that despite the fact that the quasiparticle states are a little bit 
influenced by the change of the boundary conditions, the canonical states 
remain unchanged. 

When a rather exotic nucleus is studied, it is important to check how well the 
results are converged with respect to the various approximations which are 
made in the numerical treatment of the HFB problem. First of all, we limit 
the calculation to a maximum number of partial wave given by J < Jmax, 
which has to by consistent with the mass of the nucleus. The convergence 
of the kinetic, pairing, and total energies are shown in the left part of Fig. 2. 
Between 2 Jmax = 31 and 37 the total energy difference is 1.3 keV. The changes 
in the kinetic and pairing energies are a httle bit bigger (7keV and GkeV, 
respectively); since these quantities are not at their extrema, they react more 
rapidly to a small change of the wave function. 

As discussed in section 5.3, the introduction of the cut-off radius it!cut in the 



28 



Boundary condition 




3l>1 /9 

Ir J- / ^ 


3)9q /o 

Ir Oj Z 


2 \^ /9 
J Of ^ 


y/ z 


Dirichlet 


E (MeV) 


0.926 


0.937 


1.330 


1.556 




N 


0.221 


0.577 


0.254 


0.438 




£ (MeV) 


1.169 


1.051 


1.445 


1.572 






0.251 


0.619 


0.264 


0.438 


Neumann 


E (MeV) 


0.919 


0.932 


1.328 


1.556 




N 


0.211 


0.562 


0.251 


0.438 




e (MeV) 


1.168 


1.051 


1.445 


1.572 






0.251 


0.619 


0.263 


0.438 



Table 6 

Single particle neutron levels in ^^'^Sn using the effective force Sljy4P~^^P. Calcu- 
lations have been made in a 30 fm box with an integration step of 0.2 fm. E and 
E are the quasiparticle energies and norms of the lower components while e is the 
canonical energy (54) and is the occupation factor of the corresponding canonical 
state (55). 




Fig. 2. Kinetic (top panels), pairing (middle panels) and total energies (lower panels) 

for ^^''Sn as functions of 2Jinax in a 30 fm box with an integration step of 0.2 fm 
(left) and of the pairing-field cut-off radius Rent defined in Eq. (77) (right). 

pairing channel, Eq. (77), can improve the numerical stability of the iterative 
procedure when solving the HFB equation. The right part of Figure 2 displays 
the evolution of the kinetic, pairing, and total energies with the same param- 
eters as in the left part of the Figure, except for the value of 2 J^ax which is 
fixed at 31. It can be clearly seen that beyond a rather small distance, the 
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total energy and its components are not affected by tlie introduction of this 
cut-off. Specifically, for i?cut > 17 fm the change in the total energy is less 
than O.OOlkeV. 

" 10 1 ^ 1 

> s 

i 10" 1 

3 10-1, 
10 2 , 
^10 3 1 

1 10-4: 

^ 10 S 

^ = 
■ 10-^. 

W = ; 

10 20 30 40 50 60 

Rbox [f™] 

Fig. 3. Difference between tlie total energy E of -"^^^Sn (solid line) and its average 
asymptotic value E^"-^ (see text) as a function of the box size Rhox, and the analogous 
difference for the two-neutron separation energy S2n (dashed line). 

Figure 3 shows the convergence of the total energy E and two-neutron sepa- 
ration energy as functions of the size of the box. Calculations have been 
made for -Rbox ranging from 10 fm to 60 fm with a step of 0.4 fm. The cut-off 
equivalent energy of i?cut=60MeV and the diffuseness of 1 MeV have been 
used. In order to clearly exhibit the asymptotic trend of these quantities 
we have estimated their average asymptotic values E'^"^ and ■ Since for 
-Rbox > 35 fm the results do not show any significant evolution, these val- 
ues are estimated by taking averages for 35 fm< i?box < 60 fm, which gives 
= -1 131.862 114MeV and S^^n = -1.925 548 MeV. Here, the last signif- 
icant digit corresponds to O.OOlkeV, which is the accuracy required to stop 
the iterations when we solve the HFB equations for any given value of i?box- 
The fiuctuations around the average values result from the individual contin- 
uum states entering into the cut-off window with increasing i?box- We see that 
for relatively small boxes (i? < 20 fm) the two-neutron separation energy is 
an order of magnitude more stable than the total binding energy of the nu- 
cleus. For R> 30 fm, the two quantities reach their asymptotic values with a 
random dispersion of about 0.5 keV. The value of i?box = 60 fm is obviously 
unnecessarily large but it represents a good test of the stability of the integra- 
tion procedure. From these results we conclude that the solution of the HFB 
equations with the box boundary conditions and the energy cut-off is precise 
up to about 1 keV. 

A typical asymptotic behaviour of the particle and pairing densities for neu- 
trons is shown in Fig. 4 for two different choices of the boundary condition. 
The impact of the boundary conditions only shows up within about 5 fm near 
the box edge. It is small enough to have no significant effect on the calculated 
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observables. 
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Fig. 4. Neutron particle and pairing densities in ^^'^Sn for Dirichlet (solid line) 
and Neumann (dashed line) boundary conditions. The inset represents the same 
quantities on a linear scale, in that case the impact of the boundary condition 
cannot be seen. 



Neutron particle and pairing densities for different sizes of the box are shown 
in Fig. 5. The box radius of 10 fm is obviously too small and it has only 
been included to show the effect of a very small box on the densities. For 
-Rbox > 20 fm, no differences can be seen in the linear scale. An interesting 
fact is that despite the cut that has been applied for the pairing field at 
-Rcut = 30 fm, see Sec. 5.3, no consequences can be observed on the densities 
in the asymptotic region. 
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Fig. 5. Neutron particle (left) and pairing (right) densities for different box sizes. 
In the insets the dashed lines correspond to Rhox = 10 fm while the solid lines 
correspond to the other values of -Rbox that cannot be separated in the linear scale. 



8.2 Regularization of the pairing field 



In this section, we present results obtained by using the regularization method 
proposed by A. Bulgac and Y. Yu [9]. Calculations have been made with the 
same box size and integration step as in the previous section, and the partial 
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waves have been included up to J = 43/2. The effective force used is SLy4, 
sec table 1, combined with the mixed pairing force with parameters given in 
table 2. When evaluating the densities, contributions of quasiparticlc states 
are included up to the maximum equivalent energy i?max, see Eq. (47), but 
once this maximum energy is high enough the global properties of the nucleus 
do not depend on it. 
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Fig. 6. Modulus of the difference between the total energy E and its mean value 
(solid lines) and between the average neutron gap and its mean value (dashed 
lines) as functions of the energy i^max (see text) for ^^''Sn (left) and ^^'^Sn (right). 
The mean values are calculated over the interval 60 < £^max < 80 MeV. 

It appears that the stability of results is very satisfying, even for a very exotic 
nucleus. This is shown in Figure 6, where the total energy and mean neutron 
gap are displayed as functions of ii^max for ^^°Sn and ^^°Sn. For E'max greater 
than 60 MeV, where the total energy does not show any significant evolution, 
we have evaluated its asymptotic limit E^'^\ and the analogous limit of the 
neutron mean pairing gap (A)*^"\ by averaging their respective values over the 
interval ranging from 60 to 80 MeV. The results are E^'^) = 1018.529 MeV 
and (A)(") = 1.245 MeV for ^^osn, and E^"^) = 1131.492 MeV and (A)(") = 
1.499 MeV for ^™Sn. In this interval, the energies are scattered within ±16 keV 
and the gaps within ±7 keV. Since the increase of -Bmax from 60 MeV to 
80 MeV does not change the results in a significant way, the choice of £'max = 
60 MeV has been made for the rest of this study. In principle, this value should 
be readjusted in other mass region or when using a different effective force. 

Within the cut-off prescription and regularization scheme we have calculated 
the series of even-even tin isotopes by using the SLy4(''+^^) force. The left 
part of Fig. 7 displays the binding energies per particle and the right part 
the deviation between the calculated binding energies and the experimental 
ones [25] . One can see that both methods give very similar results. The neutron 
mean gap are plotted on the left part of Fig. 8. In the two sets of calculation, 
the strengths of the pairing force have been adjusted in order to give the same 
gap in ^^''Sn. Again, we do not observe here any significant change when using 
or not the regularization scheme, although the gap is slightly reduced in heavy 
tin isotopes. Finally, the right part of Fig. 8 compares the differences between 
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Fig. 7. Left figure: Binding energy per nucleon in tin isotopes. The solid line cor- 
respond to calculations with a cut-off energy of -Ecut = 60 McV in the equivalent 
spectrum and the diamond symbols to the regularization of the energy with the 
summation of densities up to Ey^ax = 60 MeV (see text). The difference between 
the two sets of results is plotted in the inset. Right figure: Differences between the 
calculated and measured binding energies shown for those tin isotopes for which 
the data are known; solid line corresponds to the introduction of a cut-off and the 
dashed line to the regularization. 

the neutron and proton rms radii plotted with respect to that in ^^^Sn; once 
again the two methods give extremely similar results. 

8.3 Neutron drip line 

As a last example we consider the two-neutron drip line for Z — 48. The case 
of Z = 50 is not interesting for our example because there the neutron pairing 
correlations vanish at the two-neutron drip line (at least for the pairing force 
considered here), so the residual coupling between the bound and scattering 
states disappears. The two-neutron drip line is defined by the vanishing Fermi 
energy, Aat = 0. In such a system, the HFB approximation leads to a fully 
continuous quasiparticle spectrum and the minimum quasiparticle energy is 
determined by the pairing correlations only [26]. In this extreme situation 
it is important to check if the box boundary conditions have any significant 
effect on the results. To this end, we have performed calculations with the box 
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Fig. 8. Left figure: Mean neutron gap in tin isotopes, the solid line corresponds to 
the introduction of a cut-off and the dashed line to the regularization. Right figure: 
Differences between the rms radii of neutrons and protons in tin isotopes compared 
with the same quantity for ^^^Sn, the solid line corresponds to the introduction of a 
cut-off and the diamond symbol to the regularization, while the difference between 
the two sets is plotted in the inset. 



size i?box varying between 10 and 35 fm and with the regularized pairing field 
corresponding to -Bmax = 60 MeV. 
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Fig. 9. Differences between the total energy, number of neutrons, neutron rms radius, 
and neutron pairing gap, and their corresponding average asymptotic values (see 
text), calculated for the Z = 48 two-neutron drip line nucleus. 

We used the same method as in Sec. 8.1 to estimate the average asymptotic 
values (for i?box>25 fm) of the total energy 095.400 501 MeV, neutron 
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number (n)(°)=117.053 163, pairing gap ( Ajv)(") 1.150 396 MeV, and neutron 
rms radius (r jv)*-'*-'=5. 557959 fm. Differences of tliese observables relative to 
the above average values are reported in Fig. 9. It is seen that the energy is 
stable up to several keV once the radius of the box is bigger than 25 fm, and 
no further significant change is obtained by enlarging it. The mean neutron 
number is very stable too, despite the large spatial extension of the neutron 
density in such a drip line nucleus. As discussed in Ref. [26], the appearance 
of a giant neutron halo is prevented by the pairing correlations (the pairing 
anti-halo effect) and the neutron rms radius is perfectly stable with increasing 
box size. The same is also true for the neutron pairing gap. 



9 Conclusion 



The code HFBRAD (vl.OO) which we introduced in this paper allows for a 
very rapid determination of self-consistent spherical ground states of nuclei. 
It treats the pairing correlations within the HFB method and implements the 
regularization of the divergent part of the pairing energy, and hence is a perfect 
tool to study nuclei near drip lines. 

Realistic calculations must, of course, take into account the effects of the 
deformation, and the codes that are able to perform such tasks have also been 
published in Refs. [1,2]. However, deformed calculations are generally much 
more time consuming and solutions within spherical symmetry will always 
remain very useful, not only for test purposes, but also whenever results for 
many different values of input parameters are called for. This later aspect 
becomes essential when fitting of the force parameters has to be performed, 
which will constitute the main future application of the present code. 
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TEST RUN INPUT 

itinput 

force = "SLY4", mesh.points = 150, 

integ.step = 0.2, it.max = 150, 

eps.energy = l.e-9, max_delta = l.e-7, 
boundary.condition = 0, xmu = 0.65, 

bogolyubov = T, T, pairing_f orce = 3, regularization 
fenucleus neutron = 100, proton = 50, j_max =39, 25 / 
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TEST RUN OUTPUT 



********************************** PR = 8 *********************************** 

1 nucleus to compute 
Skyrme Force = SLY4 
— N = 100 ~ Z = 50 

INPUT DATA 

it.max = 150 eps energ. = l.OOOE-09 max delta = l.OOOE-07 

Boundary_condition = 

Pairing force: Volume + Surface pairing Bogolyubov transf . for both isospins 
npt = 150 integ. step = 0.200 fm Rbox = 30.00 fm r_cut = 30.00 fm 



Skyrme force parameters: 

to = -2488.913000 xO = 0.834000 tl = 486.818000 xl = -0.344000 
t2 = -546.395000 x2 = -1.000000 t3 = 13777.000000 x3 = 1.354000 
gamma = 0.166667 W = 123.000000 

tO' = -283.330000 t3' = 5312.437500 gamma' = 1.000000 
hb2/2m = 20.7355300 (with 1-body cm. correction) 



xmu = 0.65 2Jmax = 39 (neutrons), 25 (protons) 

cut.off = 60.000 MeV (with diff. 1.000 MeV) energ. step = 0.120 MeV 
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11248 


-1, 


.066778 


-17, 


,689027 


1 


.430794 


0, 


,000328 


-1131 


.8629722 


0, 


.2695E- 


-07 


56 


11161 


-1, 


.066805 


-17, 


,689058 


1 


.430822 


0, 


,000295 


-1131 


.8629981 


0, 


.2287E- 


-07 


57 


11200 


-1, 


.066827 


-17, 


,689084 


1 


.430846 


0, 


,000266 


-1131 


.8630201 


0, 


. 1940E- 


-07 


58 


11209 


-1, 


. 066846 


-17, 


,689107 


1 


.430866 


0, 


.000239 


-1131, 


. 8630389 


0, 


. 1661E- 


-07 


59 


11213 


-1, 


.066862 


-17. 


.689128 


1 


.430883 


0, 


.000215 


-1131, 


. 8630548 


0, 


. 1407E- 


-07 


60 


11203 


-1, 


.066876 


-17. 


.689145 


1 


.430898 


0, 


.000194 


-1131, 


. 8630684 


0, 


. 1202E- 


-07 


61 


11201 


-1, 


.066888 


-17. 


.689157 


1 


.430911 


0, 


.000174 


-1131, 


. 8630799 


0, 


. 1020E- 


-07 


62 


11241 


-1, 


.066898 


-17. 


.689170 


1 


.430922 


0, 


.000157 


-1131, 


. 8630898 


0, 


.8715E- 


-08 


63 


11216 


-1, 


.066907 


-17. 


.689180 


1 


.430932 


0, 


.000141 


-1131, 


. 8630982 


0, 


.7434E- 


-08 


64 


11236 


-1, 


.066914 


-17, 


.689186 


1 


.430940 


0, 


.000127 


-1131, 


.8631053 


0, 


. 6255E- 


-08 


65 


11250 


-1, 


.066920 


-17, 


,689190 


1 


.430947 


0, 


,000115 


-1131 


.8631114 


0, 


.5421E- 


-08 


66 


11299 


-1, 


.066926 


-17, 


,689203 


1 


.430953 


0, 


,000103 


-1131 


.8631166 


0, 


.4555E- 


-08 


67 


11210 


-1, 


.066930 


-17, 


,689209 


1 


.430958 


0, 


,000093 


-1131 


.8631210 


0, 


.3932E- 


-08 


68 


11196 


-1, 


.066934 


-17, 


,689211 


1 


.430962 


0, 


,000084 


-1131 


.8631248 


0, 


.3321E- 


-08 


69 


11180 


-1, 


.066937 


-17, 


,689203 


1 


.430966 


0, 


,000075 


-1131 


.8631280 


0, 


.2858E- 


-08 


70 


11202 


-1, 


. 066940 


-17, 


,689218 


1 


.430969 


0, 


,000068 


-1131 


.8631308 


0, 


.2430E- 


-08 


71 


11207 


-1, 


. 066943 


-17, 


,689218 


1 


.430972 


0, 


,000061 


-1131 


.8631331 


0, 


.2067E- 


-08 


72 


11229 


-1, 


. 066945 


-17, 


,689231 


1 


.430974 


0, 


,000055 


-1131 


.8631351 


0, 


. 1734E- 


-08 


73 


11192 


-1, 


. 066946 


-17, 


,689211 


1 


.430976 


0, 


,000049 


-1131 


.8631368 


0, 


. 1492E- 


-08 


74 


11258 


-1, 


. 066948 


-17, 


.689249 


1 


.430978 


0, 


.000044 


-1131, 


.8631382 


0, 


.1261E- 


-08 


75 


11192 


-1, 


.066949 


-17. 


.689249 


1 


.430979 


0, 


.000040 


-1131, 


.8631395 


0, 


.1116E- 


-08 



o xx±9.^ — ±.uoa3'±j — ±f.oo3^*±j j..'±ovc7(3 v.ww±w — . oooj.o;70 v.xxxo 
warning: Iterations are not progressing well and pairing seems to be small... 
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Trying to force convergence 


for it = 


= 2 














76 


11211 


-1, 


.066950 


-17 


.652243 


1 


.430981 


0, 


.000002 


-1131, 


. 8631406 


0, 


. 9561E- 


-09 


TT 
1 / 


14395 


-1, 


.066951 


-17 


.648398 


1 




u ■ 




1 lOl . 




n 

V ■ 






78 


14316 


-1, 


.066952 


-17 


.653151 


1 


.430983 


0. 


.000001 


-1131, 


.8631421 


0. 


.5991E- 


-09 


79 


14417 


-1, 


.066953 


-17 


.653151 


1 




0, 


.000001 


-1131, 


.8631428 





. 6288E- 


-09 


80 


11225 


-1, 


.066953 


-17 


.631479 


1 


.430984 


0, 


.000001 


-1131 


.8631434 


o: 


.4951E- 


-09 


81 


14323 


-1, 


.066954 


-17 


. 640449 


1 


. 43U900 


0, 


.000001 


-1131 


.8631438 


, 


.3779E- 


-09 


82 


14242 


-1, 


.066954 


-17 


. 640449 


1 


.430985 


0, 


.000001 


-1131 


.8631442 


0, 


.3551E- 


-09 


83 


11168 


-1, 


.066955 


-17 


. 640449 


1 


.430986 


0, 


.000001 


-1131 


.8631445 


0, 


.2583E- 


-09 


84 


11225 


-1, 


.066955 


-17 


.623581 


1 


.430986 


0, 


.000001 


-1131 


.8631449 


0, 


.3381E- 


-09 


85 


14356 


-1, 


.066955 


-17 


.581957 


1 


.430986 


0, 


.000001 


-1131 


.8631452 


0, 


. 2942E- 


-09 


86 


14393 


-1, 


.066955 


-17 


.581957 


1 


.430987 


0, 


.000001 


-1131 


.8631454 


0, 


. 1300E- 


-09 


87 


11269 


-1, 


.066956 


-17 


. 645585 


1 


.430987 


0, 


.000001 


-1131 


.8631455 


0, 


. 1205E- 


-09 


88 


14414 


-1, 


.066956 


-17 


.606532 


1 


.430987 


0, 


.000000 


-1131 


.8631457 


0, 


. 1428E- 


-09 


89 


14353 


-1, 


.066956 


-17 


.606532 


1 


.430987 


0, 


.000000 


-1131, 


.8631459 


0, 


. 1478E- 


-09 


iter 


Hint 




Fermi Energies 




Mean 


Delta 


Total Energy 


1 


D(E)/E 1 


90 


11162 


-1, 


.066956 


-17 


.486178 


1 


— N 

.430987 


0, 


.000000 


-1131, 


. 8631460 


0, 


. 1546E- 


-09 



Fermi Ener. = 

Mean Gaps = 

Kinetic En. = 

Pairing En. = 

Pair. Kin. En. = 



Neutrons 

■ ■ "-i'.oeegseos" 

-1.43098737 
1987.30592234 
-22.61435165 
0.00000000 



Energies (in MeV) : 

E/A = -7.545754 Variat . = 
Contributions : 

Field = -4118.250650 Spin-Or. 

(Rear. = 777.050463) 

Neutrons 

Part, numbers: 100.000000 
Radii: 5.263562 



Protons 

■ ■-i7'.486i7778" 
-0.00000038 
758.56803134 
0.00000000 
0.00000000 

0.000000 
= -67.645555 

Protons 
50.000000 
4.820502 



Total 



2745.87395368 
-22.61435165 
0.00000000 



=====> Etot = -1131.863146 

Coul. = 349.004642 
Coul.Ex. = -18.231185 
Total Charge 



5.120137 



4.886434 
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